
#include <msGeneralizedCoordinates.h>
#include <msParser.h>
#include <msZMat.h>
#include <msSystem.h>
#include <msKineticOperator.h>
#include <msLagrangian.h>
#include <muParser.h>
#include <msQuantumMotion1d.h>
#include <msLinearRigidRotor.h>
#include <msTranslation.h>

using namespace atomism;

int main(int argc, char **argv) {
 
 boost::shared_ptr<msUnitsManager> unitMng = 	msUnitsManager::New();

 unitMng->set( "Angstrom Degree amu kcal/mol ps" );

 boost::shared_ptr<msZMat> zmt =  		msZMat::New(unitMng);
 zmt->set("N N 0 1.0976");
 zmt->setId("zmat");
 
 boost::shared_ptr<msSystem> system =  		msSystem::New(unitMng);
 system->addEntity(zmt);

 boost::shared_ptr<msLinearRigidRotor> rotor = 	msLinearRigidRotor::New(unitMng);
 rotor->setMmtOfInertia(system);
 rotor->setSymmetryNumber(2);
 rotor->getParameters()->setParameter("Emax",150,*unitMng);

 boost::shared_ptr<msScalarVariable> q0 = 	msScalarVariable::New("Angstrom");
 q0->set(0.5, 0.5, 5, 0.1, 0.1);
 q0->setId("q0");

 boost::shared_ptr<msGeneralizedCoordinates> gcoors = 
						msGeneralizedCoordinates::New(unitMng);
 gcoors->addVar(q0);

 boost::shared_ptr<msParser> kinfct = 	msParser::New("Angstrom");
 kinfct->setCoordinates(gcoors);
 kinfct->setExpression("q0");
 kinfct->setId("f=q0");

 boost::shared_ptr<msKineticOperator> kinop =  	msKineticOperator::New(unitMng);
 kinop->set(gcoors,system);
 kinop->addDynamicDof(*zmt,0,kinfct);

 double De=226.0;
 boost::shared_ptr<msParser> epot = 		msParser::New("kcal/mol");
 epot->setCoordinates(gcoors);
 epot->setExpression("De * ( 1 - exp( - a*( q0 - re ) ))^2");
 epot->setId("f=Morse");
 epot->setConstante("De",De);
 epot->setConstante("a", 2.2);
 epot->setConstante("re",1.1);

 boost::shared_ptr<msQuantumMotion1d> qmotion = msQuantumMotion1d::New(unitMng);
 qmotion->setEqOfMotion(kinop,epot);
 qmotion->getParameters()->setParameter("BasisSize",701,*unitMng);
 qmotion->getParameters()->setParameter("Emax",De,*unitMng);
 qmotion->computeEigenValues();	
 qmotion->freezeKinValue( kinop->kineticFunction() );   	// kinetic fct constante : freeze it and the spectrum

 boost::shared_ptr<msTranslation> trans = msTranslation::New(unitMng);
 trans->setSystem(system);

 double Cpexp[2]; Cpexp[0] = unitMng->convert(msUnit::J_mol,29.125);   Cpexp[1] = unitMng->convert(msUnit::J_mol,37.03);  
 double Sexp[2];  Sexp[0]  = unitMng->convert(msUnit::J_mol,191.789);  Sexp[1]  = unitMng->convert(msUnit::J_mol,266.891);  
 
 int i=0;
 for( double T=300; T<3001; T+=2700 ){

      double testCv_trans = trans->Cv(T)  / (3./2. * trans->getUnits()->convert(msUnit::J_mol,csts::R) ) ; 
      double testEt_trans = trans->Eav(T) / (3./2. * trans->getUnits()->convert(msUnit::J_mol,csts::R) * T ) ; 
      double testS_trans  = trans->S(T)   / ( trans->getUnits()->convert(msUnit::J_mol,csts::R) *( log( trans->Q(T) ) + 1 +3./2. ) ); 

      fprintf(stdout,"ratio_transCv(%g)=%g\n", T, testCv_trans );
      fprintf(stdout,"ratio_transEt(%g)=%g\n", T, testEt_trans );
      fprintf(stdout,"ratio_transS(%g)=%g\n" , T, testS_trans );

      double testCv_rotor = rotor->Cv(T)  /   rotor->getUnits()->convert(msUnit::J_mol,csts::R); 
      double testEt_rotor = rotor->Eav(T) / ( rotor->getUnits()->convert(msUnit::J_mol,csts::R) * T ) ; 
      double testS_rotor  = rotor->S(T)   / ( trans->getUnits()->convert(msUnit::J_mol,csts::R) *( log( rotor->Q(T) ) + 1 ) ); 

      fprintf(stdout,"ratio_rotorCv(%g)=%g\n", T, testCv_rotor );
      fprintf(stdout,"ratio_rotorEt(%g)=%g\n", T, testEt_rotor );
      fprintf(stdout,"ratio_rotorS(%g)=%g\n" , T, testS_rotor );

      fprintf(stdout,"N2_Cp(%g)=%g (exp=%g)\n",  T, trans->Cv(T)  +  rotor->Cv(T)  +  qmotion->Cv(T)  +  unitMng->convert(msUnit::J_mol,csts::R) , Cpexp[i]  );
      fprintf(stdout,"N2_S(%g)=%g  (exp=%g)\n",  T, trans->S(T)   +  rotor->S(T)   +  qmotion->S(T)   +  rotor->getUnits()->convert(msUnit::J_mol,csts::R) ,Sexp[i] );
      i++;
      }

 }


